npj Systems Biology and Applications

AbstractsO_ST_ABSBackgroundC_ST_ABSQuantitative systems pharmacology (QSP) models provide mechanistic insight into drug response but are limited by labor-intensive, expert-driven workflows. We developed an AI-assisted QSP (AI-QSP) framework that integrates large language models (LLMs) with SBML-based modeling to enable automated reconstruction, extension, and calibration of mechanistic models. MethodsThe framework was applied to a published CAR-T QSP model. The model was reconstructed in SBML and extended via LLM-guided prompts to incorporate key resistance mechanisms: T-cell exhaustion, PD-1/PD-L1 checkpoint regulation, and tumor antigen escape. Model development followed an iterative expert-in-the-loop workflow. The resulting model (21 reactions, 9 species) was calibrated to synthetic benchmark data using 19-parameter optimization. Model credibility was assessed using ASME V&V 40 and ICH M15 principles, including global sensitivity and profile-likelihood analyses. ResultsThe calibrated model reproduced benchmark dynamics with high accuracy (mean log-RMSE = 0.132). Sensitivity analysis identified CAR-T killing and bystander cytotoxicity as dominant drivers of tumor response. Profile-likelihood analysis showed 71% of parameters were practically identifiable, with remaining parameters prioritised for future data-driven refinement. ConclusionsAI-assisted QSP modeling enables reproducible, scalable model reconstruction and evolution while maintaining mechanistic transparency and regulatory alignment. This framework provides a foundation for accelerating model-informed drug development in cell and gene therapies.

Circadian rhythms are self-sustained biological oscillations that coordinate diverse physiological processes in plants, including growth, metabolism, and environmental responses. These rhythms arise from an interconnected transcriptional translational feedback network that integrates multiple entrainment cues such as light and temperature. The plant circadian clock is organized around key regulatory loops involving CCA1, LHY, PRRs, TOC1, ELF4, LUX, and other transcriptional regulators, whose coordinated interactions ensure precise and robust oscillations. In this study, we developed an ordinary differential equation based mathematical model, building upon a previous framework to incorporate additional regulatory modules and transcriptional controls that better reflect experimentally observed behaviour. To elucidate the regulatory organization of this model, we performed a multi-layered computational analysis combining four complementary approaches: (i) period sensitivity analysis to quantify how parameter perturbations influence the systems timing, (ii) phase portrait analysis to visualize dynamic interactions among key components, (iii) knockout analysis to identify parameters essential for sustained rhythmicity, and (iv) network impact analysis using composite weighted network indices to evaluate hierarchical control across the network. Together, these analyses reveal that transcriptional repression, protein degradation, and light-regulated synthesis form the dominant control mechanisms within the circadian system. The results highlight a hierarchical and robust network structure centred on the CCA1/LHY and PRRs feedback loop, with redundant modules ensuring stability under perturbations. Thus, this model provides an improved, biologically consistent framework for dissecting the dynamic architecture of the plant circadian clock and guiding future experimental validation.

Recent studies have established that the circadian clock influences onset, progression and therapeutic outcomes in a number of chronic conditions including cardio-metabolic diseases and cancer. For the latter, they also suggest that chronotherapy offers the potential to refine current treatments and improve the development of future anti-cancer medicines. Therefore, there is a need for tools to measure the functional state of the tumoural circadian clock in patients. We have previously developed a model-led machine-learning algorithm called TimeTeller which has the potential to provide such a tool. Here we demonstrate its potential for clinical relevance by a study of 1286 breast cancer patients in which we characterise the nature of the circadian clock disruption in their tumours and demonstrate a strong nonlinear association between 10-year survival and TimeTellers tumoural clock disfunction score {Theta}. We find that good tumour clock function is antagonistic to survival.

BackgroundThe epidermal growth factor receptor (EGFR) orchestrates highly context-dependent intracellular signaling networks whose architecture varies across cell types and is frequently rewired by targeted therapeutics. Systems-level reconstruction of these networks from phosphoproteomic data remains challenging because phosphorylation measurements identify signaling nodes but do not reveal the interaction paths that propagate signals between proteins. ResultsWe developed a computational framework integrating time-resolved phosphoproteomics with graph traversal algorithms to reconstruct EGFR-initiated signaling pathways across three contexts/conditions. A sign-assignment preprocessing procedure converts quantitative phosphorylation measurements into binary activation states across time points, defining a condition-specific active node set that filters the protein-protein interaction network. Breadth-First Search combined with interaction-weighted Beam Search is then applied to the STRING interaction database (v11.5) to enumerate candidate signaling paths. Applying this framework to phosphoproteomic datasets from EGF-stimulated HeLa cells, EGF-stimulated MDA-MB-468 triple-negative breast cancer (TNBC) cells, and EGF-stimulated MDA-MB-468 cells pretreated with the SHP2 inhibitor SHP099 yielded 260 paths in HeLa cells (117 unique topologies), 293 paths in MDA-MB-468 cells (155 unique), and 292 paths under SHP2 inhibition (85 unique). HeLa cells displayed a SRC-centered architecture dominated by ERBB2 and SHC1 first-hop effectors, converging on focal adhesion, HSP90 chaperone, CRKL adaptor, and integrin signaling arms. In contrast, MDA-MB-468 cells showed a PIK3CA/PTPN11 dual-axis architecture integrating direct PI3K engagement with SHP2-mediated GRB2-IRS1-ABL1 signaling. SHP2 inhibition abolished PTPN11-mediated pathways and induced PIK3CA dominance (69.2% first-hop), accompanied by compensatory ERBB3 engagement and a computationally predicted SYK/VAV1/LCP2 node set whose biological role warrants experimental validation. ConclusionsTime-resolved phosphoproteomics-guided BFS Beam Search over STRING interaction networks captures cell-type-specific EGFR signaling architectures and drug-induced pathway rewiring. This framework provides a systematic approach for transforming phosphoproteomic measurements into mechanistically interpretable signaling hypotheses specific to the cell-type-specific contexts, directly applicable to drug resistance modeling.

Difficult-to-treat (D2T) rheumatic disease affects approximately 12% of rheumatoid arthritis patients and resists sequential biologic therapy, yet no mechanistic model explains this refractoriness as a system-level phenomenon. Here we present the 3-Axis Integrative Framework (3-AIF), a six-variable ordinary differential equation system integrating mucosal tolerance, energy-gated neuroimmune danger sensing, and integrated stress response pathways coupled through Hill-function metabolic gating. Stability analysis reveals bistable dynamics with two co-existing attractors separated by a saddle point. Bifurcation analysis demonstrates fold catastrophe with hysteresis: recovery requires greater therapeutic effort than disease prevention. Sensitivity analysis identifies three dominant parameters mapping to neuroimmune activation, energy drain, and recovery capacity. Cross-disease transcriptomic consistency analysis across six public datasets (n=310, five rheumatic diseases, four tissue types) reveals compartment-specific axis dysregulation -- circulating cells show integrated stress response activation while target tissues show pathway exhaustion -- and disease-specific axis dominance patterns consistent with model predictions.

Agent-based modeling (ABM) is a computational method for predicting the emergent outcomes of interacting, autonomous individuals in a complex system. Here, ABM is used to simulate interactions between fibroblast and myofibroblast cells during idiopathic pulmonary fibrosis (IPF) in alveolar tissue microenvironments. These microenvironments are derived from histology of a healthy human lung sample and moderate- and severe-IPF lung samples. Fibroblast differentiation, cell migration, and collagen secretion in response to the spatial distribution of the cytokine transforming growth factor-beta are captured in the ABM using NetLogo software. Results are presented from one simulated year without treatment and with mechanisms representing treatment by pirfenidone and pentoxifylline, alone and in combination. A total of 180 in silico experiments are run, analyzed, and compared in a high-throughput workflow. The effects of the initial number of fibroblasts and treatment scenarios on various metrics related to collagen accumulation and collagen invasion into alveolar regions are determined. The ABM and the analysis files are shared to facilitate model reuse. By integrating computational modeling of IPF and therapeutics, this research aims to improve understanding of fibrosis progression and assess the efficacy of novel and existing treatments targeting different mechanisms to inform decision-making for IPF treatment.

Stachys affinis (Chinese artichoke) tubers contain 50-80% stachyose by dry weight, the most concentrated dietary source of raffinose-family oligosaccharides (RFOs) known. Because humans lack sufficient -galactosidase activity, stachyose transits intact to the colon where microbial fermentation yields short-chain fatty acids (SCFAs). However, the quantitative impact of stachyose-derived SCFAs on host hepatic energy metabolism has not been systematically explored using genome-scale metabolic models. Three stachyose dose scenarios (Low/Mid/High: [~]25, 50, 100 g fresh tubers) were translated to SCFA availability vectors. Hepatic metabolic responses were simulated using Recon3D (10,600 reactions) and Human-GEM (13,417 reactions) under strict hepatocyte-like media, maximizing ATP maintenance flux (ATPM). FVA across multiple optimality thresholds (90-100%) and pFBA confirmed solution robustness. One-at-a-time sensitivity analysis characterized ATPM responses to individual parameter perturbations, and a ratio sensitivity sweep across six alternative SCFA profiles assessed dependence on assumed fermentation ratios. A targeted rescue experiment addressed model-specific propionate catabolism gaps. Both models showed dose-dependent ATPM increases (Recon3D: +71 to +286%; Human-GEM: +103 to +413% above baseline), with the 19-33% inter-model gap attributable entirely to Human-GEMs functional propionate catabolism pathway. FVA confirmed near-unique optimal solutions (ATPM ranges [~]1% at 99% optimality, widening to [~]10% at 90%). Parsimonious FBA preserved identical ATPM values while reducing total flux by [~]4-14%, confirming objective robustness. SCFA ratio sensitivity across six alternative profiles showed 27- 28% ATPM variation, indicating qualitative robustness. Butyrate yielded the highest ATP per mole ([~]22) in both models; propionate sensitivity was zero in Recon3D but [~]15.25 mmol ATPM/mmol propionate in Human-GEM. Reopening propionyl-CoA carboxylase (PPCOACm) in Recon3D under strict constraints converged ATPM to within 0.3-0.7% of Human-GEM, cross-validating both reconstructions. This reproducible dual-model pipeline identifies model-specific pathway gaps and provides cross-validated predictions to guide future experimental studies of how dietary SCFAs influence hepatic ATP metabolism.

Xylella fastidiosa is a xylem-limited phytopathogen bacterium responsible for severe diseases in numerous plant species of major agricultural importance. Despite its economic impact, its metabolism remains poorly characterized due to the bacteriums fastidious growth and the limited availability of defined culture media. In this work, we reconstructed the first pangenome-based genome-scale metabolic model for X. fastidiosa, integrating the conserved metabolic capabilities of 18 strains representing five described subspecies. The resulting core metabolic model, Xfcore, was manually curated and used to investigate the metabolic potential of the species. Model simulations predict minimal nutritional requirements that guide the formulation of defined media capable of supporting biofilm formation in vitro. Analysis of the metabolic network also suggests an undescribed metabolic pathway that enables growth on acetate as a sole carbon source. Furthermore, the model predicts that X. fastidiosa could overproduce polyamines, compounds previously associated with virulence in other phytopathogens. Experimental analyses confirm the production and secretion of polyamines in several X. fastidiosa strains, providing the first in vitro evidence of polyamine production in this pathogen. These findings suggest that polyamine biosynthesis may represent an uncharacterized virulence factor in X. fastidiosa, potentially contributing to bacterial protection against host-induced oxidative stress. Overall, the Xfcore model provides a systems-level framework to explore the metabolism of X. fastidiosa, generate testable hypotheses about its physiology and virulence, and establish a basis for future strain-specific reconstructions and host-pathogen metabolic interaction studies.

Context-specific genome-scale metabolic models are critical tools for studying cellular metabolism under dynamic conditions. However, most existing methods for deriving these models are designed for steady-state settings and may fail to preserve reactions required for transient metabolic shifts, thereby limiting their compatibility with dynamic FBA. Here, we present GeNETop, a methodology for deriving context-specific GEMs designed to preserve dynamic compatibility. GeNETop integrates flux variability analysis (FVA), network topology metrics based on the Integrated Value of Influence (IVI), and transcriptomic data to identify reactions that are both flux-flexible and structurally influential. Reactions are prioritized based on variability and maximality indices, while topology and gene expression guide further refinement, reducing dependence on fixed expression thresholds. Using batch fermentation of Saccharomyces cerevisiae as a case study, we evaluate GeNETop against established methods for context-specific metabolic reconstruction. The resulting networks remain dynamically feasible across growth phases, capture key metabolic transitions, reduce non-essential reactions, and maintain computational tractability. Overall, GeNETop enables context-specific metabolic reconstructions that are compatible with dynamic simulations while maintaining computational efficiency. By overcoming key limitations of existing approaches, the method supports a more accurate representation of time-dependent metabolic processes in biotechnology and systems biology. Author summaryCellular metabolism relies on complex networks of reactions to process nutrients, generate energy, and build essential compounds for biomass. Context-specific metabolic models aim to represent only the reactions active under a given condition, improving biological realism and reducing computational complexity in flux balance analysis simulations. However, metabolic activity adapts dynamically to changing environmental conditions, and reactions that are inactive at one stage may become essential at another. Many current reconstruction methods are designed for steady-state conditions and may exclude reactions that are required during metabolic transitions, thereby limiting their ability to describe dynamic behavior. Here, we introduce GeNETop, a novel approach that refines context-specific networks by integrating multiple layers of information. GeNETop identifies the most relevant reactions by considering their flexibility, importance within the network topology, and gene activity levels. In this way, the method generates biologically meaningful models that focus on metabolic pathways relevant under dynamic conditions. We tested GeNETop on yeast fermentation, a key process in food and biofuel production. The resulting models capture metabolic changes over time and enable stable dynamic simulations, supporting improved flux balance analysis of time-dependent metabolic processes.

Alzheimers disease (AD) is a brain disease characterized by deposition of insoluble amyloid-{beta} plaque, intraneuronal neurofibrillary tangles, and cognitive dysfunction. AD can be characterized as early-onset or late-onset based on age and genetic factors. For early-onset, these genetic factors can include amyloid precursor protein (APP), presenilin-1 (PSEN1), and presenilin-2 (PSEN2). For late-onset, these can include apolipoprotein E e4 (APOE4), and the R47H variant of triggering receptor expressed on myeloid cells 2 (TREM2). Mouse models incorporating these risk factors provide critical knowledge for studying AD pathology and preclinical studies for drug development. However, these transgenic mice depend on early-onset genetic mutations and are deficient in certain AD features that are present in late-onset. Here, we developed innovative non-linear and feature selection procedures for our cross-species translation framework, Translatable Components Regression (TransComp-R), to identify transcriptomic features in mouse models predictive of human late-onset AD pathobiology. We used the cross-species computational translatability links of TransComp-R to perform computational high-throughput drug screening and identified multiple repurposable drugs for AD treatment that targeted the sleep-wake cycle. We tested these predictions in an orthogonal, prospective cohort of human subjects treated with an orexin receptor antagonist, suvorexant. We correlated conserved protein-level biomarkers from our cross-species transcriptomics model with significant reductions in phosphorylated tau in cerebrospinal fluid collected from humans treated with suvorexant. This study demonstrates the power of computational methods like TransComp-R to enhance the utility of murine disease models for discovering new therapeutic approaches for AD. One Sentence SummaryCross-species translation modeling across different mouse models reveals sleep-relevant drug mechanisms as potentially therapeutic for Alzheimers disease.

Experimental cell systems support the development of pharmacological therapies such as glucagon-like peptide-1 receptor agonists (GLP-1RAs). However, their utility in drug discovery is limited due to study variability, which complicates formation of unified conclusions based on all available data. To address this, we conducted a comprehensive analysis of the GLP-1RA exenatide, incorporating 16 new and five pre-existing mono- or co-culture studies of human liver and pancreatic models. We employed a new pragmatic model-based approach designed to handle the common situation of heterogeneous in vitro datasets with few replicates per condition. All studies are jointly explained (disagreement<{chi}{superscript 2}-limit; 542<732), thereby providing a unified conclusion based on all studies. This work links in vitro biology to clinically relevant mechanisms, such as exenatides effect on glucose-insulin interplay, and predicts previously undescribed inter-study variabilities. Independent validation confirms predictive performance (64<83). Our new integrative approach enhances the utility of experimental cell systems in preclinical drug discovery.

Agent-based models (ABMs) are widely used to study complex multiscale biological systems, particularly in cancer research. However, their high-dimensional parameter spaces, stochasticity, and computational costs pose significant challenges for uncertainty quantification, calibration, and systematic comparison of competing mechanistic hypotheses. PhysiCell has evolved into a growing ecosystem of open-source tools supporting physics-based multicellular modeling, including model construction, visualization, and data integration. However, despite these advances, systematic support for uncertainty-aware model analysis, scalable parameter exploration, and formal calibration workflows remains limited. Here, we introduce UQ-PhysiCell, an open-source Python package that enables uncertainty quantification, calibration, and model selection for PhysiCell models using a modular and scalable workflow. UQ-PhysiCell acts as a manager of PhysiCell simulation inputs and outputs, including parameters, initial conditions, rules, and MultiCellDS-compliant objects, and provides automated orchestration of large ensembles of simulations. The framework supports multiple levels of parallelism to accelerate the analysis, including the parallel execution of independent simulations, stochastic replicates, and downstream analysis tasks. UQ-PhysiCell integrates seamlessly with established Python libraries for sensitivity analysis, optimization, Bayesian inference, and surrogate modeling, allowing users to construct customized pipelines that match their modeling goals and computational resource requirements. By decoupling model execution from statistical analysis and emphasizing extensibility and reproducibility, UQ-PhysiCell lowers the barrier to applying rigorous uncertainty-aware methodologies to agent-based modeling and supports the systematic evaluation of PhysiCell models in biological and biomedical research. Author summaryWe developed UQ-PhysiCell to address a key challenge in agent-based modeling: the systematic quantification of uncertainty in complex stochastic simulations. PhysiCell is widely used to model multicellular biological systems, particularly in cancer research; however, practical tools for uncertainty analysis, calibration, and model comparison are often developed in an ad hoc manner. This makes the results difficult to reproduce and limits the ability to rigorously evaluate competing biological hypotheses. UQ-PhysiCell provides a flexible Python framework that manages the inputs and outputs of PhysiCell simulations and enables large-scale computational analysis. We designed the software to be modular, allowing users to build their own analysis pipelines and combine different methodologies for sensitivity analysis, calibration, and model selection. Rather than enforcing a single workflow, UQ-PhysiCell supports customization to match specific scientific questions and computational requirements. To make uncertainty-aware analyses feasible for computationally intensive agent-based models, UQ-PhysiCell implements multiple parallelism strategies, enabling the concurrent execution of simulations, stochastic replicates, and downstream analyses. By promoting reproducibility, scalability, and methodological flexibility, UQ-PhysiCell helps researchers move beyond single best-fit simulations toward more reliable and interpretable computational modeling.

Interpretable representations of gene expression are used to define cellular identities and the molecular programs active within cells, two related, but distinct phenomena. In the case of microglia, a cell type with high transcriptomic, functional, and morphological heterogeneity, the predominant representation of transcriptomic data presumes the adoption of distinct molecular identities, despite a lack of easily separable transcriptional states. Here, we explore alternative transcriptomic representations by comparing two single-cell analysis methods: differential expression analysis for identities and co-expression network analysis for molecular programs. For microglia, co-expression network analysis identifies highly significant functional ontologies not resolved by differential expression analysis. The identified co-expression modules are preserved across transcriptomic datasets and suggest reducible functional programs that activate and modulate depending on context. We conclude that co-expression analysis constitutes a best practice for single cell analysis of an individual cell type and describing microglia function as concurrent molecular programs offers a more parsimonious model of microglia function.

Cells must maintain stable protein levels despite the inherently stochastic nature of gene expression, as excessive fluctuations can disrupt cellular function and impair the reliability of decision-making. Regulatory mechanisms, such as negative feedback, buffer protein fluctuations. Yet, it remains unclear how fluctuations are affected by delays that depend on a molecules specific state. Here, we develop a stochastic model in which proteins are produced in bursts as inactive molecules and pass through a series of intermediate steps before becoming active. The duration of such activation delays depends on the current level of active protein, creating a state-dependent feedback loop. Our model provides explicit analytical expressions relating the delay structure and feedback strength to the variability of active protein levels, quantified using the Fano factor, and shows that state-dependent delays can reduce fluctuations below the baseline expected from simple bursty production. Stochastic simulations confirm these predictions, and incorporating negative feedback in burst production further decreases variability while keeping system behavior predictable. These results reveal how temporal and state-dependent regulation stabilizes protein expression, offering guidance for understanding natural cellular control and designing robust synthetic gene circuits.

Fitness increase, impacted by genetic and environmental traits, requires evolutionary changes in genomes or ecological changes in habitats. Whether and how habitat reconstruction can compensate for the genetic changes remains unclear. The present study offers an experimental comparison of bacterial fitness increase via evolutionary and ecological strategies to verify that habitat reconstruction has the potential to avert genetic restriction, challenging the genetic-only view of adaptation. Six finely tuned media with different combinations and five evolved lineages with different mutations, both of which allowed equivalent increases in bacterial growth rates, were achieved through active machine learning and experimental evolution, respectively. Transcriptome changes accompanied by fitness increase were limited but divergent among the six fitted media, compared to changes that were broad but similar among the five fitted genomes. The universal transcriptome reorganization in fitted genomes and media was likely driven by surrounding metabolisms, indicating escape routes for genetic changes that increase fitness through habitat reconstruction.

Developing approaches to link structure and function is an ongoing challenge in computational and structural biology. Using a systems-level framework, we present here an analysis pipeline in a Python package, mdsa-tools, that constructs network representations of structures in a time series of trajectory frames from molecular dynamics (MD) simulations. Here, we demonstrate its use on a ribosomal subsystem. The subsystem is centered on the CAR interaction surface, a "brake pad" adjacent to the aminoacyl (A-site) decoding center that tunes protein translation rates. We leverage unsupervised learning algorithms to explore the conformational landscape of behaviors visited by two versions of the subsystem (brake-on and brake-off) that differ at the codon 3 adjacent to the A-site codon. Our network representations of MD frames embody H-bond interactions between all pairwise combinations of residues in the system. By utilizing per-frame vector representations of network edges, we can apply standard clustering and dimensionality reduction methods to explore behavioral differences between the brake-on and brake-off versions of the system. K-means clustering of frame vectors revealed a striking separation of the two system versions, consistent with principal components analysis (PCA) embeddings and Uniform Manifold Approximation and Projection (UMAP) embeddings. Dissection of K-means centroids and PCA loadings highlighted H-bond interactions between residue pairs in the ribosomes peptidyl site (P site), suggesting potential allosteric signaling across the subsystem. Author summaryWith the impressive development of computational algorithms to successfully simulate the dynamics of biological molecules over time, the exploration and incorporation of systems modes of analysis is a natural next step to begin to understand the molecular dynamics behaviors that emerge from these experiments. Following the approaches of classical molecular genetics, we used a "computational genetics" paradigm where we introduced changes (mutations) in potentially important residues, changing their identities or modifying their chemical properties, and asked how the dynamic system responded to these changes, viewing the simulations as a series of movie frames of the dynamic structure over time. Starting with network representations of each frames structure, where the nodes are residues, and the edges denote H-bond interactions between the residues, we used several unsupervised machine learning algorithms to uncover behavioral changes in the different mutated versions of the system. Applied to our ribosome neighborhood, this revealed unexpected changes in behavior at the ribosome peptidyl site (P site) in response to mutating mRNA residues on the other side of the aminoacyl site (A site) codon, suggesting long-range allosteric interactions across the neighborhood.

Understanding biological pathways requires more than static diagrams. We present WayFindR, an R package that converts pathway data from WikiPathways and KEGG into graph structures using igraph, enabling computational analysis of regulatory features such as negative feedback loops. Rooted in control theory, negative feedback is essential for system stability, yet it is often underrepresented in curated pathway data. In this study, we systematically analyzed pathway information from both databases across multiple species and found that feedback loops--particularly negative ones--are rarely captured. This gap likely reflects both biological and technical challenges. Biologically, feedback mechanisms are inherently complex and often remain uncharted due to limited experimental focus. Technically, pathway databases frequently lack standardized annotations and complete representations of regulatory interactions, especially inhibitory edges that are crucial for identifying feedback. These observations underscore the need for improved data curation and consistent annotation practices to enhance our understanding of regulatory dynamics. By bridging the gap between static pathway diagrams and dynamic systems-level insights, WayFindR enables reproducible and scalable investigation of feedback regulation in cellular networks. The WayFindR R package can be downloaded from the Comprehensive R Archive Network (CRAN) (https://cran.r-project.org/web/packages/WayFindR/index.html). The processed data along with code for download can be accessed via the GitLab repository (https://gitlab.com/krcoombes/wayfindr).

Biological rhythms are governed by intricate interactions among oscillatory subsystems, yet how they balance functional demands and energy efficiency remains unclear. We present a bimodal coupling optimization strategy where physiological systems dynamically alternate between synchronized (energy-saving) and desynchronized (function-priority) coupling modes. By employing the water-filling principle developed in communications engineering, we prove synchronized heart rate(HR)-respiration oscillations maximize energy efficiency (oxygen uptake per cardiac work). Then, system modeling confirms task/stress-induced oxygen demands enhance oxygen uptake at the cost of desynchronization and reduced efficiency. Experiments reveal a 70.36% decrease in HR-respiration synchronization during arithmetic versus relaxation, enabling 4.43% higher oxygen uptake but with 11.38% lower energy efficiency. This bimodal coupling optimization strategy is also evident in pancreatic islets, with their insulin/glucagon oscillator alternating between in-phase (energy-saving) and anti-phase (rapid glucose reduction) coupling. This framework, integrating engineering and life sciences, reveals a universal regulatory principle for biological oscillatory systems.

Tumor-induced bone disease (TIBD) arises from a complex interplay between metastatic cancer cells and the bone microenvironment, creating a self-reinforcing "vicious cycle" of bone destruction and tumor growth. Experimental evidence from our group (Buenrostro et al., Bone 113:77-88, 2018) suggests that tumor cells in the bone microenvironment early in disease rely more heavily on bone-derived growth factors, such as transforming growth factor-{beta} (TGF-{beta}), to sustain proliferation than tumor cells late in disease, which may grow independently of these factors. Here, we integrate a mechanistic, population-dynamics model of tumor-bone interactions with in vivo data to test the hypothesis that inhibiting bone resorption suppresses growth of non-adapted but not bone-adapted tumors. The model includes key regulators of TIBD, including TGF-{beta}-driven tumor proliferation, parathyroid hormone-related protein (PTHrP) secretion, and osteoblast (OB)-osteoclast (OC) coupling. Parameter calibration using data from mice injected intratibially with parental (non-adapted) and bone-adapted breast cancer cells reveals distinct parameter values for each tumor type. Bone-adapted cells exhibit a higher basal division rate and reduced sensitivity to TGF-{beta}-mediated stimulation, whereas parental-derived tumor cells depend more strongly on TGF-{beta} and secrete PTHrP at higher rates to compensate for their slower growth. Model simulations reproduce the greater bone loss observed experimentally for bone-adapted tumors and predict that, for non-adapted tumors, bone destruction results from a slower but meaningful rise in OC activity and a possible moderate decline in OBs. Simulated treatment of bone-adapted tumors with the bisphosphonate zoledronic acid stabilizes bone density but has limited or highly variable effects on tumor growth. These results suggest that OC inhibition alone may be insufficient to restrain tumor expansion once tumors have adapted to the bone microenvironment. Together, these findings support the hypothesis that tumor adaptation to the bone microenvironment governs dependence on bone-derived growth factors and response to OC-targeted therapy, underscoring the value of mechanistic modeling for elucidating tumor-bone interactions and guiding tumor-type-specific treatment strategies for TIBD.

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.